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Abstract. We study the formation of bound states and three-component bright 
vector solitons in a quasi-one-dimensional spin-orbit-coupled hyperfine spin / = 1 
Bose-Einstein condensate using numerical solution and variational approximation 
of a mean-field model. In the antiferromagnetic domain, the solutions are time- 
reversal symmetric, and the component densities have multi-peak structure. In 
the ferromagnetic domain, the solutions violate time-reversal symmetry, and the 
component densities have single-peak structure. The dynamics of the system is not 
Galelian invariant. From an analysis of Galelian invariance, we establish that the 
single-peak ferromagnetic vector solitons are true solitons and can move maintaining 
constant component densities, whereas the antiferromagnetic solitons cannot move 
with constant component densities. 
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1. Introduction 

Bright soliton is a self-reinforcing solitary wave that can traverse at a constant velocity 
withont changing its shape due to a cancellation of the non-linear and dispersive 
interactions. The various systems in which solitons have been studied include water 
waves, non-linear optics, Bose-Einstein condensates (BECs), etc. [1]. Solitons have 
been observed by manipulating the non-linear interaction near a Feshbach resonance [2] 
in a BEC of ^Li [3j and ®^Rb |1]. Solitons have also been studied in binary BECs |5]. 

In a neutral spinor BEC with a nonzero hyperfine spin /, there is no spin-orbit (SO) 
coupling between the spin of the atoms and their center-of-mass motion |S]. However, 
a synthetic SO coupling can be realized in a spinor BEC by controlling the atom- 
light interaction leading to the generation of artificial Abelian and non-Abelian gauge 
potentials coupled to the atoms |7j. Solitons have been extensively studied in spinor 
BECs without SO coupling [8]. An SO coupling with equal Rashba [9] and Dresselhaus 
uni strengths was realized experimentally by Raman dressing two atomic spin states 
with a pair of lasers m- In that study, the SO coupling between two of the three spin 
components of the / = 1 state 5Si/2 of ®^Rb — the so-called pseudospin-1/2 state — 
was considered. There are other experimental studies on SO-coupled spinor BECs m 
Solitonic structures in SO-coupled pseudospin-1/2 jlSl [H] and spin-1 BECs |15] have 
also been investigated theoretically. 

In this letter, we study two types of three-component vector solitons in an SO- 
coupled spin-1 BEC in a quasi-one-dimensional (quasi-lD) trap [T^ with muliti-peak or 
single-peak structure using a mean-field coupled Gross-Pitaevskii (GP) equation. A spin- 
1 spinor BEC is characterized by two interaction strengths, namely Cq oc (oq -I- 2a2)/3 
and C 2 oc (02 — ao)/3, where oq and 02 are s-wave scattering lengths in total spin 
/tot = 0 and 2 channels respectively im. For C 2 > 0 (antiferromagnetic) the multi¬ 
peak structure emerges, whereas for C 2 < 0 (ferromagnetic) the single-peak structure 
emerges. We use variational method to determine the bright soliton solutions for the SO- 
coupled trapless BEC in each of the two domains. The appropriate variational ansatz in 
each of the domains is constructed using the solutions of the SO-coupled single particle 
Hamiltonian. The variational analysis provides the necessary and sufficient conditions 
which Co and C 2 must satisfy to obtain a stable bright soliton. We also compare the 
variational results with the numerical solution of the GP equation. 

In Ref. |T^, only antiferromagnetic multi-peak solitons for C 2 > 0 were identified 
as the bright solitons in a three-component spin-1 SO-coupled BEG. These solitons 
are time-reversal symmetric, but are not true vector solitons as they cannot propagate 
maintaining the shape of the individual components. We demonstrate that this system 
can also support ferromagnetic single-peak solitons for C 2 < 0, provided that C 0 -I-C 2 < 0. 
These solitons break the time-reversal symmetry of the Hamiltonian. Nevertheless, they 
are shown to be true vector solitons as they can propagate with a constant velocity 
maintaining the shape of the individual components. 
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2. Spin-Orbit-coupled BEC in quasi-lD trap 


We consider an SO-coupled spinor condensate in a quasi-lD trap in which the trapping 
freqnencies along the y and ^ axes {uy and Uz) are ranch larger than that along the x 
axis {(jJx) [IS]- The single particle Hamiltonian of the condensate with equal strengths 
of Rashba [9j and Dresselhaus nm SO couplings in such a quasi-ID trap is [IH] 

( 1 ) 


Px 


^0= + 

2m 

where = —ihd/dx is the momentum operator along x axis, V{x) = mu'^x’^/2 is the 
harmonic trapping potential along x axis, and is the irreducible representation of 
the X component of the spin matrix: 

/ 0 1 0 \ 

( 2 ) 


Dx = 

72 




This SO-coupling is distinct from a previous coupling [BEn] used in the study of a 
quasi-ID BEC. 

Using the single particle model Hamiltonian ([^ and considering interactions in the 
Hartree approximation, a quasi-lD [16] spin-1 BEC can be described by the following 
set of three coupled mean-field partial differential equations for the wave-function 
components V’j 
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where Cq = 2h^{ao + 2a2)/(3m/y^), C 2 = 2 h?{a 2 — ao)/(3mly^), qq and 02 are the s-wave 


yzJ 


scattering lengths in the total spin /tot = 0 and 2 channels, respectively, pj = 
with j = 1,0,—1 are the component densities, p{x) = Yl]=-iPj ^^e total density, 
and lyz = \Jh/{mujyz) with ujyz = ydiyW is the oscillator length in the transverse y — z 
plane. For the sake of simplicity, let us transform (|^-Q into dimensionless form using 


t = uj, a; = ^, (j)j{x,t) = ^ijj{x,i), 


( 5 ) 


where Iq = \Jh/(mUx) is the oscillator length along x axis, and N is the total number 
of atoms; 
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where = ^ 2 / 2 , 7 = hkr/{mUxlo), cq = 2 A^(ao + 2 a 2 )/o/( 3 /y^), C 2 = 2 N{a 2 -a 2 )lo/{Sly^), 
Pj = |0jp with j = 1,0, —1, and p = X]j=-i The total density is now normalized 
to unity, i.e., J^^p{x)dx = 1. We present the scaled variables without tildes in 
the rest of the letter for notational simplicity. For a non-interacting trapless system 
\V{x) = Co = C 2 = 0], there are two linearly independent solutions of the SO-coupled 
set of equations @-0 with the lowest energy E m\n = —N'y‘^/2: 
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where wave functions <Fi and $2 are normalized to unity. Hence, the most general 
solution of Eqs. @-0 for a non-interacting trapless system with a hxed density n is 
given by the linear superposition of \/n^i and \/n ^2 


f \ 
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= -v/n(ai$i a2<F2), 
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where |aip -|- |a 2 p = 1 to ensure that <F is normalized to unity. 
The energy of the BEC in scaled units is 
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where F is spin density vector, whose three components Ex,Ey, and E^ are dehned as 
1 


Fx — ^ [00 (01 + <p*-l) + 0o(01 + 0-l)] ) 

= ^ ('^-1 - ■#’9 + 


( 11 ) 

( 12 ) 

(13) 

Hence, for the SO-coupled Hamiltonian with its general solution given by 0 . we get 

F. = n(|a 2 p-|«iP), (14) 

Fy = E, = 0. (15) 
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Also, the magnetization Ai = J F^dz = 0 for minimum energy solutions of the single¬ 
particle SO-coupled Hamiltonian. 

Now, let us switch on the interactions; the interaction energy per particle for the 
uniform system is [ 21 ] 



If C 2 > 0, then the BEC is in the antiferromagnetic or polar phase, and the minimum 
of Cint corresponds to |q!i| = \a 2 \ = 1/a/2 leading to |F|/n = 0. In this case, the 
wave function is time-reversal symmetric. On the other hand, for C 2 < 0, the 
BEC is in the ferromagnetic phase, and eint can be minimized if |ai| = 1, \a 2 \ = 0 or 
|ai| = 0, \a 2 \ = 1, which leads to |F|/n = 1. This corresponds to the wave functions (|^ 
apart from a multiplying phase factor. These states are degenerate, violate the time- 
reversal symmetry and are mutually connected by the time-reversal operator. These 
are the only two distinct structures which emerge as the ground states in the SO- 
coupled quasi-lD BECs. In a quasi-two-dimensional BEC with Rashba or Dresselhaus 
SO coupling, there is a circular degeneracy in the energy eigen functions of the single 
particle Hamiltonian [22]. Hence, depending upon the interaction parameters, more 
than two plane waves can also superpose resulting in different types of lattice structures 
in ground state density prohles [23] . 


3. Bright solitons 


3.1. Stationary bright solitons 

Stationary bright solitons can emerge as the ground state of a spinor BEC with attractive 
interactions |131ll5j. We use variational method to determine the bright soliton solutions 
of Q-0. As has been discussed in Sec. an SO-coupled spinor BEC can have two 
types of ground states depending upon the sign of C 2 . This necessitates the use of two 
different variational ansatz in these two domains. 

Antiferromagnetic phase (c 2 > 0); Here we consider the following variational ansatz 
to determine the shape of the soliton 

/ ± cos(7a:) \ 


$ = 


\ 


sech(cra;). 


(17) 


/ 


=F\/2i sin(7a:) 

± cos(7x) 

where a is a variational parameter and characterizes the width and the strength 
of the bright soliton. The ansatz (0 corresponds to the wave function (§ with 
Q?! = 0)2 = ±1/ a /2 multiplied by the localized spatial soliton a/2sech.{ax) instead 
of -s/n. As two solutions 0 are degenerate, and a mixing between them is allowed, the 
soliton prohle could have a multi-peak structure. Noting that in the C 2 > 0 domain, for 
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|tti| = |tt2| = l/\/2, one can have other choices for the variational ansatz like 
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$ = 
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±i sin(7a;) 


) \ 
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sech(cra;), 


(18) 


(19) 


etc, where (18) and (19) correspond to ai = —a 2 = ±1/a/2 and ai = =F*«2 = 1/v^, 
respectively, in ([^. Substituting any of these ansatz in ( |Io| , the energy of the soliton is 

+ (jco) . (20) 


E - f (-37 




( 21 ) 


The minima of this energy occurs at 
Co 

provided cq < 0. Hence, the SO-coupled spin-1 spinor BEC can support an 
antiferromagnetic bright soliton defined by (17) [or (18) or (19)] and (21), provided 
that Cq < 0 and C2 > 0. From (17) and (21) it is evident that the wavefunction of the 
bright soliton is independent of the strength of spin-exchange interactions C2. This is 
expected since for C2 > 0, there is no contribution to the energy from the C2-dependent 
term of the SO-coupled spinor BEC. 

Ferromagnetic phase (c2 < 0): Here we consider the following variational ansatz 

( 1 \ 


<hv«r = 


ae 


l'‘)X 


2^2 


V 


1 

-72 

1 


sech(cra;). 


( 22 ) 


where a is, again, a variational parameter characterizing the width and the strength 
of the bright soliton. This variational ansatz corresponds to ai = 1,^2 = 0 in (|^ 
multiplied by the localized bright soliton a/2sech.{ax) instead of \/n. In this case 

the soliton will have a single peak. Also, the ansatz like —<hvar, or ± z<hvar are equally 
reasonable choices and correspond to ai = —1, 0)2 = 0 and oi = ±f, 02 = 0, respectively. 
Substituting (p^ in (10), the energy of the soliton is 


m 


(23) 


(24) 


E = - (- 37 ^ + + aco (TC2) ■ 

0 

The minima of this energy occurs at 
^ = -^(C 0 + C2 ), 

provided Cq -|- C 2 < 0. Hence the SO-coupled spinor BEC can have a ferromagnetic 
soliton defined by (22) and (24), provided C2 < 0 and cq -|- C 2 < 0. In this case, unlike in 
the case of an antiferromagnetic soliton, the bright soliton profile is sensitive to both cq 
and C2. 
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3.2. Moving bright solitons 


If $ is static bright solitonic solution of the coupled equations (|^-Q, then the Galilean 
invariance of these equations ensures that a soliton moving with velocity v is defined as 

t) = $(a: - vt, (25) 

where a characterizes the width and the strength of the soliton. The breakdown of 
the Galilean invariance of the SO-coupled equation can be explicitly seen by using the 
transformation x' = x + vt^ t = f, where v is the velocity of the unprimed coordinate 
system with respect to primed coordinate system, then the wavefunction $ of ([^-([^ 
should transform to as 

$(x, t) = (26) 

Now, substituting (26) in ®-0 and using d/dx = d/dx' and d/dt = d/dt' + vdjdx', 
we obtain 


dt' 


2dx'^ 




(27) 


where the terms proportional to cq and C 2 have been suppressed for the sake of simplicity 
in addition to a cr-dependent additive term which does not contribute to the dynamics. 


The presence of the extra term — 7 X 3 ,^$^^(a;', t') on the right hand side of (27) shows that 


the SO-coupled Hamiltonian is no longer Galilean invariant and the SO-coupled soliton 
solution of the GP equation will depend on its velocity v. The SO-coupled equation 


(27), in the absence of trap and interactions, has the solutions $1 and <I >2 of Eq. 


with energies E = —N{'y‘^/2 — yu) and E = —iV( 7^/2 -|- yu), respectively. For u = 0, 
the two solutions ([^ were degenerate, and this degeneracy has been removed in the 
case of the SO-coupled moving solutions. In the antiferromagnetic phase, a multi-peak 
solution was possible through a mixture of two degenerate solutions (|^ for u = 0. For 
a nonzero u, the degeneracy is removed and such a mixing is not possible. This means 
that the multi-peak soliton cannot propagate with a constant velocity maintaining its 
shape and energy. For the moving multi-peak soliton profile, the variational analysis of 


Sec. 3.1 will no longer be valid. In the ferromagnetic phase, as a mixing between the two 


degenerate solutions is not allowed, one can only have a single-peak soliton which can 
propagate with a constant velocity maintaining its shape, and the variational analysis 
presented in Sec. 


3.1 remains valid. 


4. Results and conclusions 

We numerically solve the coupled equations ([^-([^ using the split-time-step Grank- 
Nicolson method [2n 125] with real- and imaginary-time propagations. The ground 
state is determined by solving ([^-([^ using imaginary-time propagation, which neither 
conserves norm nor magnetization. Both norm and magnetization can be hxed by 
transforming the wave-function components as 

0 j(a:, r - 1 - dr) = dj(l)j{x, r). 


(28) 











after each iteration in imaginary time r = —it, where dj^s with j = 1,0,—1 are the 
normalization constants. The dj^s are dehned as |20] 


do 

di 




'No + V 4 (l - M^)NiN_i + 


1 + Ni — CqA^o 


2Ni 


d-i 


'l-M-clNo 


2N_i 


(29) 

(30) 

(31) 


and here Nj = f \(j)j{x,T)\‘^dx. These normalizations ensure simultaneous conservation 
of norm and magnetization after each iteration in imaginary time. The spatial and time 
steps used in the present work are 6x = 0.05 and 6t = 0.000125, respectively. 

We consider an SO-coupled spin-1 spinor BEC of 10000 ^^Na or atoms trapped 
in a harmonic trapping potential with Ux/{2n) = 20 Hz and Uy/{2n) = bJzj (27r) = 400 
Hz. The oscillator lengths for ^^Na with these parameters are /q = 4.69 pm and 
lyz = 1.05 pm, whereas those for ®’^Rb are lo = 2.41 pm and lyz = 0.54 pm. We 
use these values of lo for writing the dimensionless GP equations ([^-([^ for the trapped 
states, whereas for solitons lo = 4.69 pm in this letter. The scattering lengths of ^^Na 
in total spin /tot = 0 and 2 channels are oq = 2.646 nm, 02 = 2.919 nm, respectively 
|26| . resulting in cq = 241.28 and C 2 = 7.76. Similarly, the scattering lengths of ®^Rb 
are oq = 5.387 nm and 02 = 5.313 nm |26], leading to cq = 885.71 and C 2 = —4.09. In 
imaginary time propagation, we use a real Gaussian function multiplied by the solution 
of the single-particle SO-coupled Hamiltonian as the initial input for the component 
wavefunctions, i.e.. 


_,2/2 / aie'T'" + 026 -'^'" \ 

^initial = -7 72 ^ 26 -'^" 


(32) 


where |ai| = \a 2 \ = l/v^ for ^^Na and |ai| = 1, \a 2 \ = 0 for ®'^Rb. Hence, by using 
different values of |ai| and |q! 2 | in (32), one can obtain different solutions corresponding 
to the same density distribution and energy. For example, the two ground state solutions 
with = 0 for ^^Na obtained by using ai = 02 = 1/v^ and ai = l/\/2 ,02 = —1/\/2 
are shown in hguresj^a) and (b), respectively. In £gures[^a) and (b), only the non-zero 
real (77) and imaginary (X) parts of the component wavefunctions are shown. In these 
two cases, wavefunctions are either purely real or imaginary and not complex. On the 
other hand, the component wavefunctions in the ground state solution for ^^Na obtained 
by using ai = 1/ a /2, 02 = i/\/2 are complex with non-zero real and imaginary parts. 
The real and imaginary parts of the component wavefunctions in this case are shown 
in figures]^ (c) and (d), respectively. The multi-peak density profile corresponding to 
these three solutions presented in hgures [^a),(b), and (c) and (d) is the same and is 
shown in Fig. [^g). The multi-peak nature of the solution in this case is consistent 
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Figure 1. (Color online) (a) and (b) show the non-zero real (7^) and imaginary 
(I) parts of wave-function components for the antiferromagnetic BEC of ^^Na with 
N = 10000, Co = 241.28, C 2 = 7.76, and 7 = 1 . For the same interaction parameters 
and SO coupling, (c) and (d) show the real and imaginary parts, respectively, of 
complex wave-function components for the antiferromagnetic BEC of ^^Na. The three 
solutions correspond to different choices of ai and 02 in (32l. (e) and (f) show the 
real and imaginary parts, respectively, of complex wave-function components for the 
ferromagnetic BEC of ®^Rb with N = 10000, cq = 885.71, C 2 = —4.09, and 7 = 1 . 

(g) shows the density distribution corresponding to (a),(b), and (c) and (d), whereas 

(h) shows the same corresponding to (e) and (f). In this and the following figures, all 
quantities are dimensionless. 
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Figure 2. (Color online) (a) Numerical and variational densities of the wave-function 
components for a bright soliton in antiferromagnetic phase with cq = —1.2, C 2 = 0.3, 
and 7 = 1. (b) The same of ferromagnetic phase with cq = —1.5, C 2 = —0.3, and 
7=1. 


with analytic results obtained in Sec. The multi-peak solution effectively leads to 
a weak phase separation between p±i and po, here weak phase separation implies that 
there are no local minima in the total density profile m This is in contrast to the 
strong phase separation possible with the model of the SO coupling discussed in Refs. 
[I^l 20 ] . where a notch appears in the total density profile at the interface separating the 
components when 7 exceeds a critical value. The solutions illustrated in figures [^a),(b), 
(c) and (d) are time-reversal symmetric. Similarly, the real and imaginary parts of the 
complex ground state solution with = 0 for ®^Rb obtained with ai = l/-\/ 2 , 0^2 = 0 
in Eq. (32) are shown in figures [^e) and (f), which lead to the single-peak density 
distribution of figure [^h). The solution presented in figures [^e), (f), and (h) violates 
time-reversal symmetry, as there are two degenerate solutions in this case connected by 
the time-reversal operation. 


p±iM - (b) 

Poix) - 




Figure 3. (Color online) Propagation dynamics of the soliton of (a) figure (a) and 
(b) figure (b) in terms of component densities p. The inset of (a) shows the variation 
in the population of the three components with time. At t = 0, the solitons located at 
cc = —10 are set into motion by multiplying the t = 0 wave function by exp(i0.2x). 
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In order to obtain the bright solitons in SO-conpled spinor BECs, we take V{x) = 0 
in ([^-Q and consider two cases: (a) Cq < 0, C 2 > 0 and (b) Cq + C 2 < 0, C 2 < 0. In 
case (a), we consider cq = —1.2,C2 = 0.3. The numerically and variationally obtained 
bright solitons, dehned by ( [l7| ) and (21) with = 0, are shown in hguresj^a). The 
multi-peak solution in this case is time-reversal symmetric. In case (b), we consider 
Co = —1.5,C2 = —0.3. The numerical and variational solutions, dehned by (22) and 
(24), in this case are shown in hgurej^b). The single-peak solution in this case breaks 
time-reversal symmetry of the Hamiltonian. It is evident from hgurej^that there is an 
excellent agreement between the numerical and variational results. 

In order to study the dynamics of the moving solitons, we hrst generate 
the stationary solitons numerically using imaginary-time propagation for both 
antiferromagnetic and ferromagnetic interactions. In order to set these solitons into 
motion with a constant velocity v = 0.2, we multiply the wavefunction components for 
the stationary soliton with exp(i0.2a:), and then use real-time propagation to study its 
evolution. We observe that in the case of the antiferromagnetic soliton, there is spin¬ 
mixing dynamics due to which the component densities are not conserved as the soliton 
moves. This is evident from hgure |^a) and its inset, which show the dynamics of the 
antiferromagnetic soliton initially located at x = —10 and the spin-mixing dynamics, 
respectively; the interaction parameters are the same as those in hgure [^a) . At f = 0 
the soliton is set into motion at a constant velocity. As the soliton moves component 
densities keep on changing without any change in the total density. On the other hand, 
if one starts with the ferromagnetic soliton at f = 0, the component densities and hence 
the total density do not change while the soliton is moving. This shown in hgure |^b) 
for the soliton initially located at x = —10 and with the same interaction parameters as 
hgure [^b). This is consistent with the analytic results of Sec. 3.2 


m 


5. Summary 

We study the generation and propagation of a vector soliton with three components in an 
SO-coupled spin-1 BEC with either antiferromagnetic or ferromagnetic interactions. In 
the antiferromagnetic case, the solutions are time-reversal symmetric and the component 
densities have multi-peak structure. In the ferromagnetic case, the solutions violate 
time-reversal symmetry and the component densities have single-peak structure. The 
GP equation for this system is not Galelian invariant. From an analysis of the Galelian 
invariance of this equation, we establish that the single-peak ferromagnetic SO-coupled 
solitons can move with constant component densities and are true solitons, whereas 
the multi-peak antiferromagnetic SO-coupled solitons change the component densities 
during motion. 
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